New Normal - Shifting Climate Baselines
Adam A. Kemberling
Updated on: 2021-05-05
Shifting Norms - Moving Climate Norms to the New Climatology Period
# load average global climatologies as annual averages
yr_avgs <- str_c(box_paths$oisst, "daily_climatologies/yearly_means/")
clim_82_avg <- stack(str_c(yr_avgs, "avg_clim_1982to2011.nc "))
clim_91_avg <- stack(str_c(yr_avgs, "avg_clim_1991to2020.nc"))
clim_shift_avg <- stack(str_c(yr_avgs, "avg_clim_shift_82to91baselines.nc"))
# # load daily_climatologies
# clim_dailys <- str_c(box_paths$oisst, "daily_climatologies/")
# clim_82_daily <- stack(str_c(clim_dailys, "daily_clims_1982to2011.nc "))
# clim_91_daily <- stack(str_c(clim_dailys, "daily_clims_1991to2020.nc"))
# clim_shift_daily <- stack(str_c(clim_dailys, "clim_shift_82to91baselines.nc"))
# climatology_lists
grid_list <- list(
"1982-2011" = clim_82_avg,
"1991-2020" = clim_91_avg,
"Climate Shift" = clim_shift_avg)
# rotate them
grid_list <- map(grid_list, raster::rotate)# convert annual_avgs to stars objects
clim_82_st <- st_as_stars(grid_list$`1982-2011`)
clim_91_st <- st_as_stars(grid_list$`1991-2020`)
clim_shift_st <- st_as_stars(grid_list$`Climate Shift`)
# Warp to world projection
wgs_84 <- st_crs(4326)
robinson_crs <- st_crs(54030)
world_moll <- st_crs("+proj=moll")
clim_82_moll <- warp_grid_projections(clim_82_st, projection_crs = "world moll")
clim_91_moll <- warp_grid_projections(clim_91_st, projection_crs = "world moll")
clim_shift_moll <- warp_grid_projections(clim_shift_st, projection_crs = "world moll")
# Get temp limits that fit both
overall_max <- max(c( max(clim_82_st[[1]], na.rm = T) , max(clim_91_st[[1]], na.rm = T) ))
overall_min <- max(c( min(clim_82_st[[1]], na.rm = T) , min(clim_91_st[[1]], na.rm = T) ))
temp_limits <- c(overall_min, overall_max)
# Map the two periods
plot_period <- function(in_stars_obj, period_label, crs_projection){
period_plot <- ggplot() +
geom_stars(data = in_stars_obj) +
geom_sf(data = world, fill = "gray90") +
geom_sf(data = st_transform(world, crs_projection), fill = "gray90") +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = temp_limits) +
guides("fill" = guide_colorbar(title = expression("Average Annual Sea Surface Temperature"~~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
coord_sf(crs = crs_projection) +
map_theme +
labs(title = period_label)
return(period_plot)
}
# Plotting the two periods
plot_82 <- plot_period(clim_82_moll, period_label = "1982-2011 Climatology", crs_projection = world_moll)
plot_91 <- plot_period(clim_91_moll, period_label = "1991-2020 Climatology", crs_projection = world_moll)
# Set diverging limits for scale
shift_limits <- max(abs(clim_shift_st[[1]]), na.rm = T) * c(-1,1)
# Plot the shift in the two
plot_shift <- ggplot() +
geom_stars(data = clim_shift_moll) +
geom_sf(data = st_transform(world, world_moll), fill = "gray95") +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = shift_limits) +
guides("fill" = guide_colorbar(title = expression("Sea Surface Temperature Difference"~~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(crs = world_moll) +
labs(title = "Global Shifts in Annual SST from Climatology Transition")
# Just plot them
stacked_clims <- (plot_82 + theme(legend.position = "none")) / plot_91
stacked_clims
plot_shift
Region Extents
The following figures look specifically at two regions. 1.) The Gulf of Maine, and 2.) the Northeastern United States’ Continental Shelf. Both of these areas have been mapped below for clarity on what data is included for any regional metric.
# load shapes/timeseries paths
gmri_areas <- get_timeseries_paths(region_group = "gmri sst focal areas")
nelme_regions <- get_timeseries_paths(region_group = "nelme regions")
# get region shapes
gom <- read_sf(gmri_areas$apershing_gulf_of_maine$shape_path)
ne_shelf <- read_sf(nelme_regions$NELME$shape_path)
# lat/lon limits
map_lims <- list(x = c(-78, -63), y = c(35, 47))
# Gulf of Maine
gom_map <- base_map +
geom_sf(data = gom,
alpha = 0.4,
color = gmri_cols("gmri blue"),
fill = gmri_cols("gmri blue")) +
scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
map_theme +
labs(title = "Gulf of Maine")
#NE Shelf
shelf_map <- base_map +
geom_sf(data = ne_shelf, alpha = 0.4,
color = gmri_cols("teal"),
fill = gmri_cols("teal")) +
scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
map_theme +
labs(title = "Northeast Shelf")
# Mapping Them together
gom_map | shelf_map
Average Regional Temperatures
# masking function
mask_sf <- function(in_ras, sf_mask){
r1 <- crop(in_ras, sf_mask)
r2 <- mask(r1, sf_mask)
return(r2)}
# Get Average Temp and Change in Climate
shift_table <- imap_dfr(grid_list, function(in_ras, climatology){
# Get means
gom_mean <- mask_sf(in_ras, gom) %>% cellStats(mean, na.rm = T)
ne_mean <- mask_sf(in_ras, ne_shelf) %>% cellStats(mean, na.rm = T)
# build table
grid_shift <- data.frame(
Climatology = rep(climatology, 2),
region = c("Gulf of Maine", "Northeast Shelf"),
avg_temp = c(gom_mean, ne_mean))
return(grid_shift) }) %>%
mutate(Climatology = factor(Climatology, levels = c("Climate Shift", "1991-2020", "1982-2011")))
# Plot the shifts
ggplot(shift_table, aes(x = avg_temp, y = Climatology)) +
geom_segment(aes(xend = 0, yend = Climatology), color = gmri_cols("gmri blue")) +
geom_point(size = 2, color = gmri_cols("orange")) +
facet_wrap(~region) +
geom_vline(xintercept = 0, color = "gray30", linetype = 3) +
labs(y = "", x = expression("Average Temperature"~~degree~C)) 
Changes to Seasonality
Loading Climatology and Heatwave Info
# Load regional timeseries
gom_ts <- read_csv(gmri_areas$apershing_gulf_of_maine$timeseries_path, col_types = cols()) %>%
mutate(time = as.Date(time))
ne_shelf_ts <- read_csv(nelme_regions$NELME$timeseries_path, col_types = cols()) %>%
mutate(time = as.Date(time))
# Side by side - Heatwave Detection
# Run climatology for both periods, merge and compare
hw_comparison <- function(region_timeseries){
# Run heatwave detection using new climate period
heatwaves_82 <- pull_heatwave_events(region_timeseries,
threshold = 90,
clim_ref_period = c("1982-01-01", "2011-12-31")) %>%
select(time, sst, seas_82 = seas, anom_82 = sst_anom, hw_thresh_82 = mhw_thresh, cs_thresh_82 = mcs_thresh,
hwe_82 = mhw_event, cse_82 = mcs_event)
# Run heatwave detection using new climate period
heatwaves_91 <- pull_heatwave_events(region_timeseries,
threshold = 90,
clim_ref_period = c("1991-01-01", "2020-12-31")) %>%
select(time, sst, seas_91 = seas, anom_91 = sst_anom, hw_thresh_91 = mhw_thresh, cs_thresh_91 = mcs_thresh,
hwe_91 = mhw_event, cse_91 = mcs_event)
# Subtract old from the new - set up flat date for plots
base_date <- as.Date("2000-01-01")
heatwaves_out <- left_join(heatwaves_82, heatwaves_91, by = c("time", "sst")) %>%
mutate(anom_shift = anom_91 - anom_82,
clim_shift = seas_91 - seas_82,
upper_shift = hw_thresh_91 - hw_thresh_82,
lower_shift = cs_thresh_91 - cs_thresh_82,
hwe_change = case_when(
hwe_82 == T & hwe_91 == F ~ "No Longer HW",
hwe_82 == F & hwe_91 == T ~ "New HW",
hwe_82 == T & hwe_91 == T ~ "Still HW",
hwe_82 == F & hwe_91 == F ~ "Still No HW",
T ~ "missing case"),
shift_direction = ifelse(clim_shift >= 0, "Temperature Increase", "Temperature Decrease"),
year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
return(heatwaves_out)
}
# Run comparisons
gom_comp <- hw_comparison(region_timeseries = gom_ts)
ne_comp <- hw_comparison(region_timeseries = ne_shelf_ts)Plotting Seasonality
# Plotting the seasonal variations
seasonality_check <- function(comp_data){
# Pull climatology
seasonal_pattern <- comp_data %>%
filter(between(time, as.Date("2019-01-01"), as.Date("2019-12-31"))) %>%
distinct(flat_date, .keep_all = T)
# Plot 82 clim
p_82 <-ggplot(seasonal_pattern, aes(x = time)) +
geom_line(aes(y = seas_82), color = gmri_cols("gmri blue")) +
geom_line(aes(y = hw_thresh_82), linetype = 2, color = "gray60") +
geom_line(aes(y = cs_thresh_82), linetype = 2, color = "gray60") +
labs(x = "", y = expression("SST"~~degree~C),
subtitle = "1982-2011 Climatology") +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
# Plot 91 Clim
p_91 <- ggplot(seasonal_pattern, aes(x = time)) +
geom_line(aes(y = seas_91), color = gmri_cols("teal")) +
geom_line(aes(y = hw_thresh_91), linetype = 2, color = "gray60") +
geom_line(aes(y = cs_thresh_91), linetype = 2, color = "gray60") +
labs(x = "", y = expression("SST"~~degree~C),
subtitle = "1991-2020 Climatology") +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
# Plot the difference
p_diff <- ggplot(seasonal_pattern, aes(x = time, y = clim_shift)) +
geom_segment(aes(yend = 0, xend = time, color = shift_direction)) +
scale_color_gmri(reverse = T) +
labs(x = "",
y = expression("Temperature Shift"~~degree~C),
color = "Shift in Climate Norm",
subtitle = "Difference in Climatologies") +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = "bottom")
return(
list(
"clim82" = p_82,
"clim91" = p_91,
"clim_shift" = p_diff)
)
}Gulf of Maine
# Return the plots
gom_seasonality <- seasonality_check(gom_comp)
# Assemble stack
(gom_seasonality[[1]] + labs(title = "Gulf of Maine")) / gom_seasonality[[2]] / gom_seasonality[[3]]
Northeast Shelf
# Return the plots
ne_seasonality <- seasonality_check(ne_comp)
# Assemble stack
(ne_seasonality[[1]] + labs(title = "Northeast Shelf")) / ne_seasonality[[2]] / ne_seasonality[[3]]
Heatwave Events
# Plotting how anomalies are now suppressed
heatwave_checks <- function(comp_data){
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(comp_data$anom_shift)) * c(-1,1)
# Base date for rectangle fill for NA
base_date <- as.Date("2000-01-01")
# Assemble heatmap plot
heatwave_heatmap <- ggplot(comp_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date,
xmax = base_date + 365,
ymin = min(comp_data$year) - .5,
ymax = max(comp_data$year) + .5,
fill = "gray75",
color = "transparent") +
# Tile for sst colors
geom_tile(aes(fill = anom_shift)) +
# points for heatwave events
geom_point(data = filter(comp_data, (hwe_82 == TRUE | hwe_91 == TRUE)),
aes(x = flat_date,
y = year,
color = hwe_change),
size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
# Colors
scale_fill_gradient2(low = gmri_cols("gmri blue"), high = "lightblue") +
scale_color_manual(values = c("black", "red")) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = "Change in Anomaly Magnitude",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "inches"),
frame.colour = "black",
ticks.colour = "black"),
color = guide_legend(title = "Impact on Heatwave Record",
title.position = "top",
title.hjust = 0.5,
override.aes = list(size = 3))) +
labs(x = "",
y = "") +
theme_classic() +
theme(legend.position = "bottom")
# Assemble pieces
return(heatwave_heatmap)
}Gulf of Maine
heatwave_checks(gom_comp) +
labs(title = "Impact of New Climatology - Gulf of Maine")
Northeast Shelf
heatwave_checks(ne_comp) +
labs(title = "Impact of New Climatology - Northeast Shelf")